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Spin-S* bilayer Heisenberg models (nearest-neighbor square 
lattice antiferromagnets in each layer, with antiferromagnetic 
interlayer couplings) are treated using dimer mean-field the- 
ory for general S and high-order expansions about the dimer 
limit for 5 = 1, 3/2, ... ,4. We suggest that the transition 
between the dimer phase at weak intraplane coupling and the 
Neel phase at strong intraplane coupling is continuous for all 
S, contrary to a recent suggestion based on Schwinger boson 
mean-field theory. We also present results for S — 1 layers 
based on expansions about the Ising limit: In every respect 
the S = 1 bilayers appear to behave like 5* = 1/2 bilayers, 
further supporting our picture for the nature of the order- 
disorder phase transition. 

PACS numbers: 



I. INTRODUCTION AND MEAN-FIELD 
ARGUMENTS 

The bilayer Heisenberg antiferromagnet, 
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where the intralayer couplings run over nearest neigh- 
bors on square lattices and the S are S = 1/2 quan- 
tum spins, has turned out to be an excellent testing 
ground for notions of quantum criticality. A variety of 
controlled calculations, such as strong-coupling expan- 
sions^at zero temperature (starting with the work._pf 
HidaEl and extended by two of the present authorsa 
finite-temperature quantum Monte-,Carlo calculations^ 
and high-temperature expansionsQ fit extremely well 
into the general field-theoretic framework for zero- 
temperature order-disorder critical points in unfrustrated 
two-dimensional quantum antiferromagnets. 

The critical point in the S = 1/2 bilayer Heisen- 
berg model is between a Neel ordered phase when A = 
Ji / J 2 is sufficiently large (greater than the critical value 
A c = 0.394) and a magnetically disordered "dimer" phase 
which is the adiabatic continuation of the A = ground 
state (which is a product of interlayer singlet pair wave- 
functions). The universality class of the phase transition 
is the same as for finite-temperature transitions in d = 3, 
0(3) symmetric models. 

It seemed most natural to us that for larger values 
of the quantum spin S the same scenario would play 



out: The value of A c would change, but qualitatively 
nothing would be different. It ^therefore came as a sur- 
prise when Ng, Zhang, and Mac! described the results of 
Schwinger boson mean-field (SBMF) calculations, which 
indicate that for S > S c the transition between magnet- 
ically ordered and disordered phases is first-order. The 
calculation leads to S c « 0.35, and hence indicates that 
for all physical values of S the transition should be first- 
order; but it is reasonable to accept the suggestion of 
Ng et al. that SBMF theory leads to phase diagram of S 
versus A which is qualitatively valid even if S c is under- 
estimated. Two more salient points regarding the SBMF 
results should be noted. At large S, the value of A at the 
dimer-Neel phase boundary scales as 1/5; and at suffi- 
ciently large S the ground state on the dimer side of the 
boundary is the pure interlayer singlet pair state (i.e., 
identical to the ground state at A = 0). 

Here we will argue that the results of SBMF theory 
for phase transitions in the bilayer Heisenberg model are 
qualitatively incorrect, even in the limit of large S where 
one might expect such an approximation to be reliable. 

Ng et al. offer an elementary argument in support 
of their proposed scenario, which it is informative to 
reconsider. They note that the ideal Neel state and 
the interlayer singlet state become degenerate for A = 
Ai = 0(1/5) (with the interlayer singlet state energet- 
ically favored for A < Ai), while linear spin- wave the- 
ory indicates that the magnetization should vanish at 
A = A 2 = 0(1/S 2 ) (with Neel order existing for A > A 2 ). 
At large S, it is clear that Ai > A 2 . Ng et al. suggest 
that, on decreasing A from infinity, the transition that 
one might suspect should take place at A 2 (presumably 
a continuous transition from the Neel to dimer phase) 
is pre-empted by a first-order transition at Ai. This is 
certainly consistent with their SBMF calculations. 

However, we suggest an alternative interpretation. We 
believe that the argument in favor of a first-order tran- 
sition at Ai is too simplistic: Using simple trial states 
to represent the ground states in the respective phases is 
dangerous, particularly since it is a comparison of sub- 
leading (in powers of S) terms in the variational en- 
ergies for each phase which leads to the estimate for 
Ai. Furthermore, there is an elementary calculation 
which strongly suggests that the dimer phase is unsta- 
ble with respect to the development of Neel order for 
A > A 3 = 0(1/S 2 ). Consider a single dimer, subject 
to a staggered magnetic field which is produced self- 
consistently by the staggered magnetization of the sys- 
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tern. Let J2Xo(S) denote the staggered susceptibility 
of a single dimer. Then, within this "dimer mean- field 
approximation" one has A3 = l/ixo- One can calcu- 
late xo directly from second-order perturbation theory, 
based on explicit formulae for 3-j symbols,cl to find that 
A3 = 3/165(5 + 1). As promised, the domain of stabil- 
ity of the dimer phase is 0(l/5 2 ), in agreement with the 
spin-wave estimate for the stability of the Neel phase. 

It should be clear that the dimer mean-field approxi- 
mation underestimates the stability of the dimer phase, 
just as the standard Weiss theory overestimates Curie 
temperatures. One could reasonably object to our ar- 
gument, then, by suggesting that the corrections to this 
mean-field theory could be large: in particular, A3 might 
be pushed up to a value of order 1/5. In the follow- 
ing section, we describe high-order perturbation expan- 
sions about A = for bilayers with 5=1, 3/2, . . . , 4. 
Combined with the existing results for 5 = 1/2, they 
strongly suggest that the corrections to dimer mean-field 
theory are order unity rather than order 5. In addition, 
for 5=1 bilayers we present a wide variety of results 
based on expansions about A = and about Ising mod- 
els. They are all qualitatively similar to what is found 
for 5=1/2 bilayers, providing additional confirmation 
that the phase diagram is unchanged with increasing 5 
(in particular, there is no evidence for any phase besides 
dimer or Neel) and offering accurate results which could 
prove useful in interpreting data from experimental real- 
izations of 5 = 1 bilayers. 

II. SERIES EXPANSIONS AND 
EXTRAPOLATIONS 

A. Triplet Gap and Antiferromagnetic Susceptibility 

Perturbation expansions!^ for the ground-state energy 
Eq, the triplet excitation spectrum, and the antiferro- 
magnetic susceptibility x have been carried out for bi- 
layers with all integer and half-integer values of 5 from 1 
to 4. The series coefficients for integer 5 bilayers for the 
minimum gap A, corresponding to wave vector (n,ir), 
and x are presented in Table Q. For 5 = 1/2 bilayers 
such expansions have already been presented by Zhengtl 
toO(A n ). 

Fig. [l] is a set of "scatter plots" of estimated criti- 
cal A values versus estimated critical exponents derived 
from unbiased Dlog-Pade approximants to the gap se- 
ries for 5 = 1, 2, 3 and 4 bilayers. An analogous plot 
for 5 = 1/2 is presented in Ref. ||. In every case they 
have very similar character, with nearly all exponent es- 
timates lying slightly above the anticipated value of 0.71, 
and with clear correlations between critical point and 
exponent values. (Scatter plots for the antiferromagnetic 
susceptibility have very similar character; likewise for the 
gap and susceptibility of half-integer 5 bilayers.) As 5 
increases the series actually become better behaved, in 



that the approximants derived from a fixed number of 
terms are more tightly clustered (neglecting the outliers) 
and estimated exponents are closer to the expected value. 
Biasing the exponent one can obtain more precise esti- 
mates of the critical points A c , as listed in Table |n[ The 
trend is clear: with increasing 5, the critical value of A 
is approaching a constant multiple of the dimer mean- 
field value 3/165(5+ 1). If we plot 1/A C as a function 
of R = 5(5 + 1) (shown in Fig. ||), we can see that the 
results can be remarkably well fitted by a straight line 

1/A C = c(R - R c ) (2) 

with an intercept at 1/A C = slightly above R = 0. The 
constants c and R c determined by a linear least squares 
fit are: 

c = 3.72(1), R c = 0.068(2) . (3) 

For general 5 we have found empirically, based on the 
ground state energy Eq, the minimum triplet gap A, and 
the antiferromagnetic susceptibility x series for the eight 
values of 5 that have been explicitly calculated, that the 
coefficients of A" can be expressed as polynomials of order 
n in the variable R. For example, to third order 

A = 1 - (8/3)RX + [(8/9)i? - (32/27)i? 2 ]A 2 (4) 
+[(56/135)i?- (116/135)i? 2 - (3824/1215)i? 3 ]A 3 . 

Further coefficients are listed in the Appendix. One can 
then define r = XR and consider the R — ► 00, r — ► const 
limit. The resulting series in r corresponds to the terms 
in the double-series (in R and A) of the form R n X n . Its 
Dlog-Pade approximants are as well behaved as the series 
for the larger values of 5 shown above, with unbiased ap- 
proximants clustering and exponent estimates just above 
the anticipated value. Biasing the critical point estimates 
leads to a critical r of 0.2691(7) (this is consistent with 
the value of c in Eq. @)), so the ratio of the exact critical 
point value in the large-5 limit to that coming from the 
dimer mean-field theory is 1.435(4). 

We can also try to use the double-expansions in R and 
A to analytically continue to unphysically small values 
of 5. An interesting feature of the general-5 double- 
series (both for the A and x) is that all terms of the 
form R°X n - 1 appear to have vanishing coefficients. Con- 
sequently, for 5 = there is no phase transition to a 
Neel-ordered state. Presumably there is a critical spin 
Sci less than 1/2, at and below which the dimer state is 
stable for all values of A. It has not been possible to ob- 
tain a reliable estimate of 5 C by constructing Dlog-Pade 
approximants at fixed values of 5: for 5 < 0.35 the ap- 
proximants depend strongly on the number of terms used. 
However, the extrapolation from larger 5, summarized in 
Eqs. (§) and (§), suggests that 5 C = 0.064(2). 
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B. Further results for S = 1 bilayers 

Although we have presented evidence above that the 
dimer phase is unstable for A larger than a critical value 
of 0(1/ 5 2 ), we did not firmly establish that the struc- 
ture of the phase diagram is simply dimer phase-critical 
point-Neel phase. The fact that the unbiased exponent 
values are consistent with a Lorentz-invariant, d = 3, 
0(3) universality class is certainly suggestive. For the 
S = 1 bilayers we have done much more. In addition 
to expansions about A = we have constructed expan- 
sions about Ising models, generating a set of results anal- 
ogous to those presented by Zhengq for 5=1/2 bilayers. 
Taken together, these leave little room for doubt about 
the phase diagram. 

Because the calculations follow those presented in 
Rcf. U so closely we refer the reader to Sec. Ill of that 
paper for a description of the calculation (note that the 
parameter y in that paper corresponds to 1/A). Expan- 
sions were obtained to 0(x 9 ) (with x denoting the ratio 
of transverse to longitudinal exchange strength) for the 
sublattice magnetization M, uniform transverse suscep- 
tibility X-Li and triplet excitation spectrum. From these 
we could derive the gap to the optical branch of the spin- 
wave spectrum at wave vectors (0,0) (A opt , the mini- 
mum gap) and (0,7r) (Ax), spin-wave velocity v, and 
spin- wave stiffness p s . The expansion coefficients will 
not be presented here but are available from the authors 
on request. Results of series extrapolations are shown 
in Figs. [| f|, and ||. They demonstrate that expansions 
about the Neel state lead to estimates of the domain of 
stability of the Neel phase which are completely consis- 
tent with the critical point found by expanding about 
A = 0. This is nicely confirmed by Fig. |6[ which displays 
the ground-state energy per site. The values from the 
dimer expansion match very smoothly on to those from 
the Ising expansion around the critical point, whereas for 
a first-order transition there would be a discontinuity in 
slope at the transition. All of the results presented in this 
subsection represent the best available theoretical values 
for experimentally accessible properties of Neel-ordered 
5 = 1 bilayers, which should prove useful in the inter- 
pretation of experimental data if any such systems are 
studied. (In real compounds single-ion anisotropy is al- 
ways present to some degree. The present calculations 
could readily be extended to include such terms in the 
Hamiltonian.) 

Finally, it is amusing to note that in the dimer phase 
for 5 > 1/2 bilayers, sufficiently close to A = 0, there ex- 
ist spin-2 elementary excitations. Let us consider 5=1 
bilayers specifically. Triplet spectra for various values of 
A in the dimer phase are shown in Fig. ^ while quintuplet 
spectra are shown in Fig. ||. (These spectra are obtained 
by direct sums of the terms in the series to the maxi- 
mum order available.) The latter spectra lose physical 
significance when the quintuplet excitations can decay 
into multiple triplet excitations. One might imagine this 



happens for arbitrarily small values of A, since at A = 
the quintuplet gap is three times the triplet gap. How- 
ever, the parity, with respect to interchange of layers in 
the bilayer system, of the quintuplet excitations is op- 
posite to that of the triplet excitations, and symmetry 
forbids decay of the quintuplet excitations into an odd 
number of triplet excitations. Hence stable spin-2 exci- 
tations lie between the 2-triplet and 4-triplet continua for 
sufficiently small A. 

III. DISCUSSION 

The order-disorder transition in bilayer Heisenberg 
antiferromagnets appears to be a problem for which 
Schwinger boson mean-field theory leads to qualitatively 
incorrect results even in the limit of large 5. We have 
shown by high order perturbation expansions about the 
limit of uncoupled interlayer singlets that there is strong 
evidence for continuous transitions between dimer and 
Neel phases with critical values of J1/J2 scaling as 
1/5(5 + 1). There seems to be no-,evidence in favor of 
the scenario described by Ng et aZ.,cl in which the transi- 
tion would become first-order for sufficiently large 5. Of 
course this does not imply that the phase diagram for 
bilayer Heisenberg antiferromagnets is completely uni- 
versal. If further-neighbor interactions or higher-order 
exchange interactions (for example, terms in the Hamil- 
tonian of the form (Sj • Sj) 2 ) were allowed then a wide 
variety of new phases could be stabilized. However, for 
the simplest bilayer models it appears that simple argu- 
ments based on the instability of the dimer phase to Neel 
order (using dimer mean-field theory) and the instability 
of the Neel phase to a spin-disordered phase (using linear 
spin- wave theory to determine where the staggered mag- 
netization vanishes), which both suggest a critical point 
at J\/J% = 0(l/5 2 ), are correct. 
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APPENDIX: 

The double expansions in A and R = 5(5+1) for 
the ground-state energy per site Eq/N, triplet minimum 
energy gap A, and antiferromagnetic susceptibility x f° r 
spin- 5 bilayers are presented in full here. 

E _ -R 2\ 2 R 2 X 3 R 2 
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+A J R 0.166019988 - 4.50447076 R 



+8.81922354 R 2 + 5.13123542 R 3 - 14.6722347 R 4 

+X 6 R ^0.112022751 - 5.43708747 R + 24.41683462 R 

-23.43145101 R 3 + 17.32806918 R 4 - 21.97846687 R 5 

+X 7 R (0.091548205 - 6.26174454 R + 45.07198387 R 2 
-73.28433315 R 3 + 29.56605425 R 4 + 59.68596039 R 5 
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FIG. 1. Plots of critical A values versus critical exponents 
for Dlog-Pade approximants to the triplet gap series. All 
approximants that use all available terms in the series (circles) 
and all but the last term (squares) are displayed. Note that 
many of the approximants are difficult to distinguish on these 
plots, but the "outliers" represent single approximants. 
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(A2) FIG. 2. Plot of the inverse of critical point (1/A C ) versus 
S(S + 1): the crosses with error bars are the estimates bi- 
ased by the critical index, and the dashed line is a linear fit: 
1/Ac = 3.72{S(S + 1) - 0.068]. The insert enlarges the region 
near S = 0. 



FIG. 3. Rescaled energy gaps A opt and Ax as a function 
of Jil(J\ + J 2)- The solid curves at large Ji/{J\ + J2) are 
extrapolations based on dimer expansions; the crosses with 
error bars are the estimates from Ising expansions; the long 
dashed lines at small Ji/{J\ + J2) are the results of the linear 
spin-wave theory;til and the open circles connected by a sho+4 
dashed line are results from the theory of Millis and Monienli3 
(with p s , an input to their theory, taken from Fig. ^). The 
vertical dotted line indicates the critical point. 
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FIG. 4. The staggered magnetization M and uniform 
(A3) perpendicular susceptibility xx f° r S = X bilayers versus 
(./2/Ji) 1 ^ 2 as estimated by Ising expansions. 



FIG. 5. The spin-wave velocity v and the spin-wave stiff- 
ness p s versus (J2/Ji) 1//2 as estimated by Ising expansions. 
At the critical ratio J2/J1 ~ 7.18, the estimate of v from 
dimer expansions is also shown. 
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FIG. 6. The rescaled ground-state energy per site 
Eq/(AJ\ + J2)/N as a function of J2HJ1 + ^2)- The crosses 
with dashed line connecting them are the results from the ex- 
pansion about the Ising limit, while the solid line gives the 
results from the dimer expansion. 

FIG. 7. Plot of the spin-triplet excitation spectrum 
A(k x ,k y ) along high-symmetry cuts through the Brillouin 
zone for the system with coupling ratios J1/J2 = 0.05, 0.1, 
0.12, 0.1393 (shown in the figure from the top to the bottom 
at (n,n) respectively). The lines are the estimates by direct 
sum to the dimer series, and the points (circles with error bar 
for the case of J1/J2 = 0.1393 only) are the estimates of the 
Pade approximants to the dimer series. 

FIG. 8. Plot of the quintuplet excitation spectrum 
A q (k x ,k y ) along high-symmetry cuts through the Brillouin 
zone for the system with coupling ratios J1/J2 = 0.05, 0.1, 
0.12, 0.1393. The lines are the estimates by direct sum to 
the dimer series, and the points (circles with error bar for the 
case of J1/J2 = 0.1393 only) are the estimates of the Pade 
approximants to the dimer series. 
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TABLE I. Series coefficients c n for dimer expansions of the minimum tripiet gap A = J2 Y2 n c ™-^ 
and the antiferromagnet susceptibility x = J^ 1 S c «-^ n f° r ^ — 1> 2, 3, and 4 bilayers. 





S = 1 


S = 2 


S = 3 


S = 4 






Minimum triplet gap 







1.0000000000 


1.0000000000 


1.0000000000 


1.0000000000 


1 


-5.3333333333 


— 1.6000000000X 10 1 


— 3.2000000000X 10 1 


-5.3333333333X10 1 


2 


-2.9629629630 


—3.7333333333 x 10 1 


— 1.6000000000X 10 2 


-4.5629629630X10 2 


3 


-2. 7786008230 xlO 1 


-7.0826666667X 10 2 


— 5.5573333333X 10 3 


-2.5514008230X10 4 


4 


-4. 0140832190 xlO 1 


—3.5921481481 x 10 3 


— 6.1199238095X 10 4 


-4.8582313672 x10 s 


5 


-3.3454379922 xlO 2 


- 1.0569742866 x 10 s 


— 3.5299272466X 10 6 


-4.6061398131 xlO 7 


6 


-1.0532000998X10 3 


— 9.1597242219X 10 s 


—6. 1760036073 XlO 7 


-1.3547279280X10 9 


7 


-9.7019942999X10 3 


-2.6882225528X 10 7 


— 3.6317201466X 10 9 




8 


-3.0018452174X10 4 








9 


-3.0081157280x10 s 












Antiferromagnet ic 


susceptibility 







2.6666666667 


8.0000000000 


1.6000000000 xlO 1 


2.6666666667X10 1 


1 


2. 8444444444 xlO 1 


2.5600000000X10 2 


1.0240000000 xlO 3 


2. 8444444444 xlO 3 


2 


2.4800000000 xlO 2 


6.7120000000X10 3 


5.3728000000 xlO 4 


2.4880000000 x10 s 


3 


2.0721234568X10 3 


1.7017800000x10 s 


2. 7320720000 xlO 6 


2.1109167901 xlO 7 


4 


1.6499774192 xlO 4 


4.1278260778X10 6 


1.3301975846x10 s 


1.7154007839X10 9 


5 


1.2891136900x10 s 


9.8615644145X10 7 


6.3850802835 xlO 9 


1.3748488989X10 11 


6 


9.9036806734 x10 s 


2.3183667236x10° 


3.0164597503X10 11 


1.0845538535 xlO 13 


7 


7. 5397998426 xlO 6 


5.4074404771 xlO 10 






8 


5.6895542718X10 7 









TABLE II. The critical point A c for S — 1/2 to S = 4 estimated by biasing the critical index 
(where the case of S = 1/2 is already reported in Ref. ^|), and its ratio to that predicted by dimer 
mean-field theory. The results of a linear fit, = 3.72[S(S + 1) - 0.068], are also listed. 
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1/2 
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3/2 


2 


5/2 


3 


7/2 


4 


Ac 


0.3942(8) 


0.1393(7) 


0.07303(10) 


0.0453(1) 


0.03098(7) 


0.02254(5) 


0.01716(5) 


0.01349(5) 


A fit 


0.3942 


0.1391 


0.07301 


0.04532 


0.03096 


0.02253 


0.01714 


0.01349 


A c /A 3 


1.58 


1.49 


1.46 


1.45 


1.44 


1.44 


1.44 


1.44 
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10 12 
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0.8 
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x- 
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0.08 - 
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0.06 
0.04 
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(J 2 /Jl) 
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